% =================================================================
% 主程序：使用LU分解求解二维变系数泊松方程
% 方程形式：-∇·(a(x,y)∇u) = f(x,y)
% 功能说明：
% 1. 生成二维泊松方程刚度矩阵
% 2. 构建空间变系数矩阵a(x,y)和高斯源项f(x,y)
% 3. 使用带主元选择的LU分解法求解线性系统
% 4. 可视化数值解、系数分布和源项分布
% =================================================================

% 关闭所有图形窗口
close all

% ========================= 输入参数设置 ==========================
n = 20;         % 单方向内部节点数
a_amp = 12;     % 系数函数a的幅度参数
f_amp = 1;      % 源项f的幅度参数
x_0 = 0.5;      % 高斯分布中心x坐标
y_0 = 0.5;      % 高斯分布中心y坐标
c_x = 1;        % x方向标准差系数
c_y = 1;        % y方向标准差系数
h = 1/(n+1);    % 网格步长计算（区域[0,1]×[0,1]）

% ====================== 矩阵与向量生成 ==========================
% 生成n²×n²的二维泊松方程刚度矩阵（五点差分格式）
S = DiscretePoisson2D(n);

% LU分解（带部分主元选择）
[L, U, P] = lu(S);

% 构建系数矩阵a(x,y)（高斯型变系数）
C = zeros(n,n);
for i = 1:n
    for j = 1:n
        C(i,j) = 1 + a_amp*exp(-((i*h-x_0)^2/(2*c_x^2) + ...
                  (j*h-y_0)^2/(2*c_y^2)));
    end
end

% 创建对角矩阵D（用于系数处理）
D = zeros(n^2,n^2);
for i = 1:n
    for j = 1:n
        D(j+n*(i-1), j+n*(i-1)) = C(i,j);
    end
end

% 构建源项向量f（高斯分布源）
f = zeros(n^2,1);
for i = 1:n
    for j = 1:n
        f(n*(i-1)+j) = f_amp*exp(-((i*h-x_0)^2/(2*c_x^2) + ...
                       (j*h-y_0)^2/(2*c_y^2)));
    end
end

% ====================== 方程求解过程 ==========================
% 计算右端项b = D^{-1}f
b = zeros(n^2,1);
for i = 1:n
    for j = 1:n
        b(n*(i-1)+j) = f(n*(i-1)+j)/C(i,j);
    end
end

% 前向替换求解 Lv = P*b
% v = ForwSub(L, P*b);
v = L\(P*b);
% 后向替换求解 Uw = v
% w = BackSub(U, v);
w = L'\v;
% 最终解计算 u = h²*w
u = h^2 * w;

% ====================== 结果可视化 ==========================
% 将解向量转换为网格格式（含边界零点）
Z = zeros(n+2,n+2);
for i = 1:n
    for j = 1:n
        Z(i+1,j+1) = u(j+n*(i-1));
    end
end

% 生成网格坐标
x1 = 0:h:1;
y1 = 0:h:1;

% 解场二维等高线图
figure(1)
surf(x1, y1, Z)
view(2)
colorbar
xlabel('x_1'), ylabel('x_2'), zlabel('u(x_1,x_2)')
title(['数值解场分布 系数幅度=', num2str(a_amp), ', 网格数=', num2str(n)])

% 解场三维表面图
figure(2)
surf(x1, y1, Z)
colorbar
xlabel('x_1'), ylabel('x_2'), zlabel('u(x_1,x_2)')
title(['三维解场视图 系数幅度=', num2str(a_amp), ', 网格数=', num2str(n)])

% 系数a(x,y)分布图
Z_a = zeros(n+2);
for i = 1:n+2
    for j = 1:n+2
        Z_a(i,j) = 1 + a_amp*exp(-(((i-1)*h-x_0)^2/(2*c_x^2) + ...
                    ((j-1)*h-y_0)^2/(2*c_y^2)));
    end
end
figure(3)
surf(x1, y1, Z_a)
xlabel('x_1'), ylabel('x_2'), zlabel('a(x_1,x_2)')
title(['系数分布 a(x,y) 幅度=', num2str(a_amp)])

% 源项f(x,y)分布图
Z_f = zeros(n+2);
for i = 1:n+2
    for j = 1:n+2
        Z_f(i,j) = f_amp*exp(-((x1(i)-x_0)^2/(2*c_x^2) + ...
                    (y1(j)-y_0)^2/(2*c_y^2)));
    end
end
figure(4)
surf(x1, y1, Z_f)
xlabel('x_1'), ylabel('x_2'), zlabel('f(x_1,x_2)')
title(['源项分布 f(x,y) 幅度=', num2str(f_amp)])
